Liquid-liquid equilibrium for monodisperse spherical particles 
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A system of identical particles interacting through an isotropic potential that allows for two preferred 
interparticle distances is numerically studied. When the parameters of the interaction potential are 
adequately chosen, the system exhibits coexistence between two different liquid phases (in addition 
to the usual liquid-gas coexistence). It is shown that this coexistence can occur at equilibrium, 
namely, in the region where the liquid is thermodynamically stable. 
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I. INTRODUCTION 

Traditionally, liquid-liquid separation for single com- 
ponent systems was not considered as a real possibility, 
because of the rapidly changing liquid structure, which 
smoothly varies on temperature and pressurea. This view 
is changing. Experimentally, a whole family of network 
forming fluids, including some tetrahedrally coordinated 
materials, is suspected to have transitions between liquid 
statcsH. Usually this transition occurs well inside the su- 
percooled region, where the liquid is already in_fl_glassy 
state, and then it is related to polyamorphismEm. The 
numerical modeling of these materials shows that inter- 
action potentials can be constructed for which there are 
transitions between different liquid states. For instance, 
the potentials used to mimic the properties of water gen- 
erate two liquid phases in the supercooled regiona. The 
main characteristic of these materials (and of the interac- 
tion potentials used to model them) which is responsible 
for their anomalous properties is the possibility of gen- 
erating open, low coordinated network structures at low 
pressure, while collapsing to more compact structures at 
higher pressuresa. 

This key ingredient is captured by a class of model 
spherical potentials, the so-called soft core Stell-Hemmcr 
potentialscl. They include a strict hard core at some dis- 
tance ro and a soft core at a larger distance r\. In this 
way, two typical distances between particles occur, and 
there is a collapse from the largest to the smallest dis- 
tance on applying pressure. The properties of these sys- 
tems give a remarkably consistent description of most 
of the anomalous behavior of tetrahedrally coordinated 
materialsHcl. 

The possibility of liquid-liquid separation has been ex- 
tended in the last time by the finding that it can oc- 
cuc. even at thermodynamic equilibrium. Katayama et 
al.u report that phosphorus (which is a network form- 
ing fluid) shows liquid-liquid coexistence at equilibrium, 
providing the first experimental evideace of this phe- 
nomenon. In addition, Glosli and Ree£2l find numerical 
evidence of liquid- liquid coexistence at equilibrium in car- 
bon. This was the first direct numerical evidence of the 



phenomenon for a single component system. The inter- 
particle potential used to model carbon-carbon interac- 
tion is non-isotropic, and incorporates information of the 
internal (electronic) degrees of freedom of the atoms. It 
is shown that one of the two coexisting liquids is mainly 
sp coordinated, whereas the other shows a dominant sp 3 
character. 

On the basis of these results, it is natural to ask 
whether some kind of soft-core potentials (with no in- 
ternal degrees of freedom for the particles) can produce 
an equilibrium liquid-liquid transition. We will show, by 
providing concrete examples, that soft core Stell-Hemmcr 
potentials are able to do that, then showing that the phe- 
nomenon may occur also for systems that, in addition to 
being monodisperse, have particles which are spherical 
and 'rigid', i.e., in cases where internal degrees of free- 
dom play no role. This is of basic interest since this model 
is much simpler, and then more transparent to analyze 
than real examples of material having liquid-liquid co- 
existence at equilibrium. It also provides a benchmark 
for the study of tetrahedrally coordinated materials with 
liquid-liquid critical points in the supercooled region. 
The relation between the existence of the liquid-liquid 
critical point and the anomalous thermodynamic proper- 
ties of these materials can be studied here in a context 
of thermodynamic equilibrium. On a different context, 
our results make plausible the existence of liquid-liquid 
separation at equilibrium in particular cases of colloidal 
dispersions, which is of potential technological interest. 



II. THE MODEL 

We study a model of particles interacting through a 
potential which is the sum of attractive and repulsive 
parts V(r) = Va{t) + Vr{t). The repulsive part is given 
as 



V R (r) 

V R (r) = e (ri - r) / (n - r ) 
V r {t) 



00 for r < tq 
for tq < r < n 
for r > n, (1) 
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namely, a hard core at distance ro, and a linear (repul- 
sive) ramp between ro and n. The attractive part is 
taken also to be linear in r, of the form 

V A {r) = -7 (r 2 - r) / (r 2 - r ) for r < r 2 

FA(r) = for r > r 2 (2) 

(r 2 > ri, 7 > 0). The repulsive part of the potential 
favors one of the two distances ro or r\ between neighbor 
particles depending on the value of external pressure and 
the attraction intensity 7. Interparticle distance t\(tq) is 
favored for pressures lower (higher) than some crossover 
value Pq. The crossover pressure (for 7 = 0) is given 
roughly as Pq — Eq/(vi — vq), where vq and v\ are the 
specific volumes of the states with nearest particles at 
distances ~ ro or ~ ri . For P near Pq the system has an 
anomalously large compressibility. The attractive part 
of the potential may turn this anomaly into a first order 
transition through a van der Waals mechanismB. 

We will show results for two-dimensional (2D) and 
three-dimensional (3D) systems. The parameters we use 
are: ri = 1.72r , r 2 = 4.8r , and 7 = 0.27e for 2D and 
?"i = 1.72ro, r 2 = 3.0ro, and 7 = 0.31eo for 3D. Temper- 
ature will be measured in units of £o/fcs, and pressure in 
units of £0/^0 m 2D, and £0/^0 m 3D. Calculations were 
done by standard Monte Carlo techniques, in systems 
with periodic boundary contitions. Some of the results 
to be presented where done in the NVT ensemble, where 
the volume of the system is kept fixed at each simulated 
value, and pressure is calculated through standard for- 
mulae. Other results correspond to the NPT ensemble, 
where external pressure is fixed, and the volume of the 
system is taken as an additional variable during the sim- 
ulation. 



III. RESULTS 

A. Two dimensions 

Fig. [l| shows the isotherms we get from the simula- 
tions for a system of 200 particles in the 2D case, near 
the position in which we find a liquid-liquid critical point 
for our parameters, namely Tll — 0.10, Pll ~ 0.07, 
and vll ~ 1.85^. We see in fact that when reducing 
T below ~ 0.1, the isotherms get a loop, typical of a 
first order transition. This transition separates two liq- 
uid phases that we call (following the notation in water) 
low density liquid (LDL) and high density liquid (HDL). 
We have checked that the position and characteristics of 
this transition do not change when going to a system of 
1000 particles. As an example, we show in Fig. ^ the 
T = 0.105 and T = 0.0875 isotherms of systems of 200 
and 1000 particles. Isotherms are seen to coincide with 
the 200 particles case, except within the coexistence re- 
gion. Here the isotherm is more flat in the case of larger 



system. This is in fact what expected, since in an in- 
finitely large system isotherms are strictly flat in the co- 
existence region. Fig. || shows snapshots of the system 
at T = 0.0875, within the LDL and HDL phases and in 
the coexistence region. The HDL (LDL) phase is char- 
acterized by a larger amount of particles at distance ro 
(ri) from their neighbors. The snapshot at coexistence 
[Fig. ||(b)] shows indeed that different neighborhoods of 
high density and low density coexist in the system. The 
existence of a loop in the isotherms below the critical 
temperature, instead of a flat region, even when the sys- 
tem clearly separates in different regions corresponding 
to the two different phases [Fig. ||(b)] is due to the non- 
vanishing contribution to the total energy of the surface 
energy between the two phases. This contribution only 
vanishes in infinite systems. 

From the isotherms of Fig. [j] we see that the LDL has 
a density anomaly, since there is a temperature range 
in which volume increases as temperature is reduced, at 
constant pressure. This density anomaly (which is rem- 
iniscent x>f the similar phenomenon in water) is known 
to occurEl in core-softened models even in the absence of 
an attractive part of the interparticle potential, namely, 
when U a = 0. 

We will now present the equilibrium P-T phase di- 
agram of the system, close to the liquid-liquid critical 
point, to show clearly that the liquid-liquid coexistence 
we are observing occurs in thermodynamic equilibrium, 
and it is not preempted by the appearance of any other 
phase (gaseous or crystalline). 

The P-T phase diagram is presented in Fig. and 
the different borders between phases were obtained as 
explained in the following. The position of the equilib- 
rium first order line between LDL and HDL phases in the 
P-T phase diagram was determined from the isotherms 
of Fig. [| and from the values of the enthalpy h = e + Pv 
(not shown) by thermodynamic integration, calculating 
differences in Gibbs free energy G from the relation 

^-gi= f^dP-^dT (3) 
T 2 T\ J 1 T T 2 y ' 

where 1 and 2 stand for two set of values Pi , T\ and P 2 , 
T 2 and the integration is through an arbitrary (reversible) 
path in th P-T plane. In this way we can compare G on 
both sides of the first order line, performing an integra- 
tion surrounding the critical point, and the first order line 
can be determined with good precision. This procedure 
is more reliable than the standard Maxwell construction, 
although the results we have obtained in this case are 
almost equivalent. 

Concerning the crystalline structures, we know in fact 
that the present model shows a variety of different crys- 
talline structures, . wh ich are stable in different regions 
of the P-T planeBlii When testing for the stability of 
these structures, the fact that some of them take a very 
long time to appear spontaneously in simulations that 
start in the fluid phase is usually a problem, and then the 
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equilibrium coexistence line between these phases and the 
fluid has to be determined in some other way. We then 
applied the procedure already used in [11] consisting in 
the inclusion of an additional external potential, with an 
spatial periodicity corresponding to the crystalline lattice 
to be studied. The existence of this term (with a strength 
given by some parameter W) allows to construct a path 
in the P-T-W space that takes the system smoothly (i.e., 
without crossing any phase transition) from the fluid to 
the crystalline phase, and the difference in free energy 
AG between the two phases at some fixed values P = Pq 
and T = Tq can be calculated by a generalization of ex- 
pression (Q), namely 

A(G/T)\ Po , To =J^ d P + ^ d W-^dT, (4) 

where ew is the potential energy per particle in the ex- 
ternal potential. In order to move the system from the 
crystalline to the fluid phase, the integration consists typ- 
ically in the following close path: increasing W from zero 
to some large value; increasing T from T to a large value; 
decreasing W down to 0; decreasing T down to To. In 
this way, the difference in free energy between the fluid 
phase and the possible crystalline phases was determined 
at fixed points in the P-T phase diagram. Then AG can 
be calculated for the whole plane by thermodynamic inte- 
gration of the pure phases. Using this procedure we have 
verified that the only crystalline phase that is thermodi- 
namically stable within the P-T sector shown in Fig. [| is 
the triangular structure, with a lattice parameter ~ ro. 
Other crystalline structures appear only for lower tem- 
peratures of higher pressures. 

Note in Fig. ^ that the melting of the crystal is anoma- 
lous in the region near liquid-liquid coexistence, i.e, both 
liquid phases are denser than the solid around this point, 
which is indicated by the negative slope of the liquid- 
crystal coexistence line. This is probably not a general 
rule, but is the case we found in all simulations we per- 
formed. In this figure the line of maximum density is 
also indicated. Note that this line does not end exactly 
at the critical point, but somewhere on the first order 
line, within the LDL region. 

At very low pressure the condensed (liquid or solid) 
phases are unstable against the formation of a gaseous 
phase, as in a standard fluid. The condensed-gas co- 
existence line can be numerically determined in a reli- 
able way by the standard simulation methods we use 
only close to the (liquid-gas) critical point, where the 
metastability range is small. Typical isotherms around 
this point are shown in Fig. ^, from which the loca- 
tion of the liquid-gas critical point is estimated to be at 
P LG ~ 0.02, T LG ~ 0.37. Using the position of this crit- 
ical point, we estimate the liquid-gas coexistence line by 
fitting to a van der Waals equation, chosen to give the 
correct values of Tlc and Plg- This is the line that is 
seen in the right part of Fig. [| (note the change in the 
horizontal scale in this sector). At temperatures (~ 0.1) 



where we observe liquid-liquid coexistence, the transition 
pressure to the gas phase is of the order of 10~ 6 . Then 
we are sure that the liquid-liquid transition occurs in a 
region in which liquid is thcrmodynamically stable with 
respect to gas. 

The crystal-liquid- liquid triple point we obtained turns 
out into a gas-liquid-liquid triple point (still at equilib- 
rium) by increasing the intensity of the attraction be- 
tween particles. In fact, this moves all coexistence lines 
to lower pressures compared to the condensed-gas line, 
and at a certain point the crystal-liquid-liquid triple point 
is forced to be in the metastable region with respect to 
the gas phase. At this point, the gas-condensed line, 
together with the liquid-liquid line determine a liquid- 
liquid-gas triple point. For instance, increasing 7 from 
0.27 to 0.285 produces a liquid-liquid-gas triple point lo- 
cated atP~10~ 6 ,T~0.10. The liquid-liquid critical 
point is in this case located at P ~ 0.03, T ~ 0.12. 
At still higher values of the attraction the whole liquid- 
liquid line moves into the metastable P < region, and 
the equilibrium liquid- liquid transition is lost. 

B. Three dimensions 

For the 3D case, we show in Fig. ^ the isotherms ob- 
tained for a system of 300 particles. They again show the 
existence of the liquid- liquid critical point. The curves 
are smooth for T > 0.08, whereas there is a jump in v 
at the coexistence pressure for T lower than this value. 
Also apparent from these curves is the existence (as in 
2D) of density anomalies for the low density liquid. 

The structure factor S(k) of the 3D liquid is shown for 
the two coexisting liquids at T = 0.075 in Fig. 0. The two 
main peaks corresponding to interparticle distance ro and 
n are indicated. There is an abrupt change in the relative 
weight of these two peaks at the liquid-liquid transition, 
corresponding to an abrupt change in the mean distance 
between first neighbor particles. 

For the simulations shown in Fig. |^ we did not find 
evidence of the appearance of any crystalline structure, 
suggesting that the liquid-liquid coexistence occurs at 
equilibrium. To be sure, we should investigate the sta- 
bility of the different possible crystalline structures, as 
we did in the 2D case. For the 3D system however, this 
is a much harder work, as we actually do not know all 
the crystalline structures for this case (see [12]). Then 
we take another route, seeking for other indicators of 
stability of the liquid phase, to rule out the possibil- 
ity that we are observing a liquid-liquid transition in a 
metastable region. So we will lopk for the fulfillment 
of the Hansen- VerletEJ and LowerO criteria for crystal- 
lization. Although they stand only as approximate, and 
mainly heuristic, they have been verified in a variety 
of cases. For simple (i.e. Lennard Jones) systems, the 
Hansen- Verlet criterion states that freezing occurs when 
the first peak of S(k) reaches a (quasi-)univcrsal value 
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of rsj 2.85. A straightforward generalization of this crite- 
rion to 'anomalous' fluids tells that this relation should 
be checked for all peaks in S(kp. The results for S(k) 
of Fig. ^ show that all peaks are lower than this value, 
suggesting liquid stability. The Lowen criteria for freez- 
ing indicates that this occurs when the relation between 
long and short times single particle diffusion coefficient 
drops below a value of ~ 0.1. This dynamical criterion 
is independent of the nature of the crystalline structure 
the system freezes to. Although it has been proposed 
and verified in systems with Brownian dynamics, we have 
verified in simple cases (i.e. hard spheres) that it is also 
fulfilled with the Monte Carlo dynamics, and so we apply 
it tentatively here. The relation between long and short 
time values of the diffusivity D/Dq (calculated with the 
Monte Carlo dynamics) is shown in Fig. || along the 
T = 0.075 and T = 0.09 isotherms. The relation D/D 
is in the whole range well above the value of 0.1 expected 
at freezing. Then all indicators coincide in telling that 
the liquid-liquid coexistence we arc observing occurs in 
truly thermodynamic equilibrium. From Fig. || we also 
see that the evolution of D on density shows a maxi- 
mum which is weiLknown to occur in tetrahedrally coor- 
dinates materialaLj and it is also typical of core-softened 
potentialsta. 



IV. CONCLUSIONS 

In conclusion, we have shown that a system of identi- 
cal, spherical, and rigid particles is able to display liquid- 
liquid coexistence at equilibrium (both in two and three 
dimensions). We showed explicitly this for an interac- 
tion potential that favors two different equilibrium dis- 
tances between particles. The study of the present (or 
some related) model displaying liquid-liquid coexistence 
at equilibrium is likely to contribute to the understand- 
ing of the same phenomenon in cases where it occurs in 
the supercooled regionfl, where it is much harder to study 
due to metastability against crystallization and long re- 
laxation times. Our findings also open the possibility of 
observing liquid-liquid coexistence in monodisperse sys- 
tems of rigid and spherical particles, i.e., colloidal sys- 
tems, where interaction potentials with— more than one 
equilibrium distance are easily obtainecO. 
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FIG. 1. Isotherms from the simulation of a 2D system of 
200 particles, at temperatures as indicated (symbols are the 
simulated points, lines are guides to the eye). There is a crit- 
ical point at T ~ 0.10, P ~ 0, 07, below which the isotherms 
get a typical van der Waals loop. Note the expansion upon 
decreasing T at constant P (a density anomaly) in the results 
for w/rg>2.1. 
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FIG. 2. Comparison of isotherm obtained in simulations of 
systems with 200 and 1000 particles. They are coincident, ex- 
cept within the first order loop, where the one corresponding 
to the larger system is more flat. This is the right tendency 
since the isotherm at coexistence is absolutely flat in the limit 
of systems of infinite size. 



FIG. 3. Snapshots of the system of 1000 particles at the 
points indicated in Fig. |^. Black dots represent the strict 
hard core of the particles at ro. In (a) and (c) the system is 
almost completely in one of the two liquid phases (HDL and 
LDL, respectively). In (b) there is a clear coexistence of the 
two phases. 
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FIG. 4. Pressure-temperature phase diagram of the 
two-dimensional model showing the existence of the liq- 
uid-liquid and liquid-gas critical points. The temperature of 
maximum density (TMD) is indicated. 
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FIG. 5. Isotherms near the liquid-gas transition for a 
two-dimensional system of 500 particles. 
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